//    Copyright (C) Kherty.  All rights reserved.
using System;
using OpenLS.Spreadsheet.Extensibility;

namespace OpenLS.Spreadsheet.StandardAddIn.Statistics
{
    internal class TTest
    {
        public double PairedTTest(double[] sample1, double[] sample2)
        {
            double meanDifference = Statistician.MeanDifference(sample1, sample2);
            return
                tTest(meanDifference, 0, Statistician.VarianceDifference(sample1, sample2, meanDifference),
                      sample1.Length);
        }


        public double HomoScedasticTTest(double[] sample1, double[] sample2)
        {
            if ((sample1 == null) || (sample2 == null ||
                                      Math.Min(sample1.Length, sample2.Length) < 2))
            {
                throw new ArgumentException("insufficient data");
            }
            return HomoScedasticTTest(Statistician.Mean(sample1),
                                      Statistician.Mean(sample2), Statistician.Variance(sample1),
                                      Statistician.Variance(sample2), sample1.Length,
                                      sample2.Length);
        }

        public double HeteroScedasticTTest(double[] sample1, double[] sample2)
        {
            if ((sample1 == null) || (sample2 == null ||
                                      Math.Min(sample1.Length, sample2.Length) < 2))
                throw new ArgumentException("insufficient data");
            return HeteroScedasticTTest(Statistician.Mean(sample1),
                                        Statistician.Mean(sample2), Statistician.Variance(sample1),
                                        Statistician.Variance(sample2), sample1.Length,
                                        sample2.Length);
        }

        private static double t(double m, double mu, double v, double n)
        {
            if (v == 0)
                throw new ErrorException(Errors.DivError);
            if (n == 0)
                throw new ErrorException(Errors.DivError);
            return (m - mu)/Math.Sqrt(v/n);
        }

        private static double homoScedasticT(double m1, double m2, double v1,
                                      double v2, double n1, double n2)
        {
            double pooledVariance = ((n1 - 1)*v1 + (n2 - 1)*v2)/(n1 + n2 - 2);
            return (m1 - m2)/Math.Sqrt(pooledVariance*(1d/n1 + 1d/n2));
        }


        private static double heteroScedasticT(double m1, double m2, double v1, double v2, double n1, double n2)
        {
            return (m1 - m2)/Math.Sqrt(v1/n1 + v2/n2);
        }

        private static double tTest(double m, double mu, double v, double n)
        {
            double tt = Math.Abs(t(m, mu, v, n));
            TDistribution tDistribution = new TDistribution(n - 1);
            return 1.0 - tDistribution.CumulativeProbability(-tt, tt);
        }

        private static double HomoScedasticTTest(double m1, double m2, double v1,
                                          double v2, double n1, double n2)
        {
            double tt = Math.Abs(homoScedasticT(m1, m2, v1, v2, n1, n2));
            double degreesOfFreedom;
            degreesOfFreedom = n1 + n2 - 2;
            TDistribution tDistribution = new TDistribution(degreesOfFreedom);
            return 1.0 - tDistribution.CumulativeProbability(-tt, tt);
        }

        private static double HeteroScedasticTTest(double m1, double m2, double v1,
                                            double v2, double n1, double n2)
        {
            double tt = Math.Abs(heteroScedasticT(m1, m2, v1, v2, n1, n2));
            double degreesOfFreedomTop = v1/n1 + v2/n2;
            degreesOfFreedomTop = degreesOfFreedomTop*degreesOfFreedomTop;
            double degreesOfFreedomBottom = (v1/n1)*(v1/n1)/(n1 - 1) + (v2/n2)*(v2/n2)/(n2 - 1);
            double degreesOfFreedom = degreesOfFreedomTop/degreesOfFreedomBottom;
            TDistribution tDistribution = new TDistribution(degreesOfFreedom);
            return 1.0 - tDistribution.CumulativeProbability(-tt, tt);
        }
    }
}